
% note as currently written, the result WIDTH is put out in units of pixels
% and NOT in microns

function [r,vesselTotalIntensity,surroundInten,maxInten,validtraj,localgcamplevel] = extractWidefieldTrajParams_v1_1(trackingfn,mouseID,sessionID,dosave,gcamppixeladjustment,prestimframes,stimframes,fps)
% takes output of highResWidefieldWrapper analysis and organizes results
% neatly into a table with relevant parameters
%
% inputs:
% trackingfn - the filename of the tracking file outputted by
% highResWidefieldWrapper
%
% mouseID - character array identifying the mouse ID that that trackign
% comes from
%
% sessionID - number indicating which tracking session that data comes from
% gcamppixeladjustment - undoes the pixel correction for 4x4 vs. no binning
% binning; not relevant to any experiments after 1/23/2024

if nargin < 4 || isempty(dosave)
    dosave = false;
end
if nargin <6 || isempty(fps)
    fps = 1/.14;
end
if nargin < 7 || isempty(prestimframes)
    prestimframes = 2:10;
end
%% load in trajectories
load(trackingfn,'traj','stimID_idx','vesselID','uid','localgcampdistance','localgcamplevel')
if nargin < 5 || isempty(gcamppixeladjustment)
    temp = max(localgcamplevel,[],3,'omitnan');
    temp = max(temp,[],2);
    if prctile(temp,95) < 9.5
        gcamppixeladjustment = true;
    else
        gcamppixeladjustment = false;
    end
end
if nargin < 7 || isempty(stimframes)
    stimframes = 11:39;
end
if nargin < 8 || isempty(fps)
    fps = 1/.14;
end
numTrackedSegs = size(traj(1).vesselCent,1);
numframes = size(traj(1).vesselTotalInten,2);
numtrials = length(traj);


%% determine valid trajectories ahead:
validtraj = filterNoisyWidefieldTrajectories_v1_0(traj,prestimframes);
validtraj2 = filterWidefieldTrajByWidth(trackingfn,prestimframes);
validtraj = validtraj & validtraj2;
%% smooth tracking
vesselTotalIntensity = nan(numTrackedSegs,numframes,numtrials);
surroundInten = nan(numTrackedSegs,numframes,numtrials);
maxInten = nan(numTrackedSegs,numframes,numtrials);
vesselWidth = nan(numTrackedSegs,numframes,numtrials);

for n = 1:numTrackedSegs
    for n2 = 1:numtrials
        cur = traj(n2).vesselTotalInten(n,:);
        missing = isnan(cur);
        if sum(missing) > 0 && sum(missing)<(.25*length(cur))
            cur = resample(cur,1:length(cur));
            vesselTotalIntensity(n,:,n2) = wdenoise_NVC_v2(cur);
        elseif sum(missing) == 0
            vesselTotalIntensity(n,:,n2) = wdenoise_NVC_v2(cur);
        end
        if validtraj(n2,n)
            traj(n2).vesselTotalInten(n,:)= vesselTotalIntensity(n,:,n2);
        else
            traj(n2).vesselTotalInten(n,:)=nan;
        end

        cur = traj(n2).surroundInten(n,:);
        missing = isnan(cur);
        if sum(missing) > 0 && sum(missing)<(.25*length(cur))
            cur = resample(cur,1:length(cur));
            surroundInten(n,:,n2) = wdenoise_NVC_v2(cur);
        elseif sum(missing) == 0
            surroundInten(n,:,n2) = wdenoise_NVC_v2(cur);
        end
        if validtraj(n2,n)
            traj(n2).surroundInten(n,:)= surroundInten(n,:,n2);
        else
            traj(n2).surroundInten(n,:)=nan;
        end

        cur = traj(n2).maxInten(n,:);
        missing = isnan(cur);
        if sum(missing) > 0 && sum(missing)<(.25*length(cur))
            cur = resample(cur,1:length(cur));
            maxInten(n,:,n2) = wdenoise_NVC_v2(cur);
        elseif sum(missing) == 0
            maxInten(n,:,n2) = wdenoise_NVC_v2(cur);
        end
        if validtraj(n2,n)
            traj(n2).maxInten(n,:)= maxInten(n,:,n2);
        else
            traj(n2).maxInten(n,:)=nan;
        end

        cur = traj(n2).vesselWidth(n,:);
        missing = isnan(cur);
        if sum(missing) > 0 && sum(missing)<(.25*length(cur))
            cur = resample(cur,1:length(cur));
            vesselWidth(n,:,n2) = wdenoise_NVC_v2(cur);
        elseif sum(missing) == 0
            vesselWidth(n,:,n2) = wdenoise_NVC_v2(cur);
        end
        if validtraj(n2,n)
            traj(n2).vesselWidth(n,:)= vesselWidth(n,:,n2);
        else
            traj(n2).vesselWidth(n,:)=nan;
        end
    end
end
disp('smoothing done')

%% extract relevant parameters from gcamp:
% I think in practice we'll have to get the median response for each stim
% type in each session and use that? otherwise may be too variable...
analysisrange = 12:39;  %window of frames to query here
baselinerange = 2:5;        %because of smoothing need to take a smaller baseline range for this
distance2gcamp = median(localgcampdistance(:,analysisrange,:),2,'omitnan');
distance2gcamp = reshape(distance2gcamp,numTrackedSegs,numtrials);

totalgcamp = median(localgcamplevel(:,analysisrange,:),2,'omitnan');
totalgcamp = reshape(totalgcamp,numTrackedSegs,numtrials);
if max(totalgcamp(:)<8)
    gcamppixeladjustment = true;
end
if gcamppixeladjustment
    totalgcamp = log(exp(totalgcamp)*16);
end

prestimgcamp = median(localgcamplevel(:,baselinerange,:),2,'omitnan');
prestimgcamp = reshape(prestimgcamp,numTrackedSegs,numtrials);
if gcamppixeladjustment
    prestimgcamp = log(exp(prestimgcamp)*16);
end
%% get baseline values:
[width_baseline,vesselTtotalIntensity_baseline,surroundInten_baseline,maxInten_baseline] = widefieldTrajectoryNormalization_2led_v1_0(localgcamplevel,traj);
temp = mean(validtraj,1);
temp = temp > .75;          %trajectories where more than this threshold are valid

vesselTtotalIntensity_baseline(~temp) = nan;
surroundInten_baseline(~temp) = nan;
maxInten_baseline(~temp) = nan;
width_baseline(~temp) = nan;

%% normalize everything to its baseline:
maxInten = maxInten./repmat(maxInten_baseline,1,numframes,numtrials);
maxInten = (maxInten-1)*100;
surroundInten = surroundInten./repmat(surroundInten_baseline,1,numframes,numtrials);
surroundInten = (surroundInten-1)*100;
vesselTotalIntensity = vesselTotalIntensity./repmat(vesselTtotalIntensity_baseline,1,numframes,numtrials);
vesselTotalIntensity = (vesselTotalIntensity-1)*100;

%% add in additional smoothing step because wavelet filtering sometimes misses some higher frequency noise:
maxdeviation = 5;       %can have consistent deviation threshold since everything is normalzied
for n = 1:size(maxInten,1)
    for n2 = 1:size(maxInten,3)
        cur = maxInten(n,:,n2);
        if sum(~isnan(cur)) < 5
            continue;
        end
        cursm = smooth(cur,'sgolay',1);
        d = abs(cur-cursm');
        cur(d>maxdeviation) = cursm(d>maxdeviation);
        maxInten(n,:,n2) = cur;

        cur = surroundInten(n,:,n2);
        cursm = smooth(cur,'sgolay',1);
        d = abs(cur-cursm');
        cur(d>maxdeviation) = cursm(d>maxdeviation);
        surroundInten(n,:,n2) = cur;

        cur = vesselTotalIntensity(n,:,n2);
        cursm = smooth(cur,'sgolay',1);
        d = abs(cur-cursm');
        cur(d>maxdeviation) = cursm(d>maxdeviation);
        vesselTotalIntensity(n,:,n2) = cur;
    end
end

%% extract relevant parameters from the ios tracking
evokedMaxVessel = nan(numTrackedSegs,numtrials);
evokedMaxSurround = nan(numTrackedSegs,numtrials);
evokedMaxTotalVessel = nan(numTrackedSegs,numtrials);
evokedIntegratedMaxVessel = nan(numTrackedSegs,numtrials);
evokedIntegratedSurround = nan(numTrackedSegs,numtrials);
evokedIntegratedTotalVessel = nan(numTrackedSegs,numtrials);
time2peak = nan(numTrackedSegs,numtrials);
for n = 1:numTrackedSegs
    for n2 = 1:numtrials
        if validtraj(n2,n)      %only bother tracking trajectories that are deemed valid
            temp = extractTrajParams_v3_4(maxInten(n,:,n2),[stimframes(1),stimframes(end)],fps,false,true);
            evokedMaxVessel(n,n2) = temp.maxDilation(1);
            evokedIntegratedMaxVessel(n,n2) = temp.totalDilation;
            time2peak(n,n2) = temp.maxDilation(2);

            temp = extractTrajParams_v3_4(surroundInten(n,:,n2),[stimframes(1),stimframes(end)],fps,false,true);
            evokedMaxSurround(n,n2) = temp.maxDilation(1);
            evokedIntegratedSurround(n,n2) = temp.totalDilation;
    
            temp = extractTrajParams_v3_4(vesselTotalIntensity(n,:,n2),[stimframes(1),stimframes(end)],fps,false,true);
            evokedMaxTotalVessel(n,n2) = temp.maxDilation(1);
            evokedIntegratedTotalVessel(n,n2) = temp.totalDilation;
            
        end
    end
end





%% create a table with all these values in it
evokedMaxVessel = evokedMaxVessel(:);
evokedIntegratedMaxVessel = evokedIntegratedMaxVessel(:);
evokedMaxSurround = evokedMaxSurround(:);
evokedIntegratedSurround = evokedIntegratedSurround(:);
evokedMaxTotalVessel = evokedMaxTotalVessel(:);
evokedIntegratedTotalVessel = evokedIntegratedTotalVessel(:);
time2peak = time2peak(:);
distance2gcamp = distance2gcamp(:);
totalgcamp = totalgcamp(:);
prestimgcamp = prestimgcamp(:);
vesselID = repmat(vesselID,numtrials,1);
vesselTtotalIntensity_baseline = repmat(vesselTtotalIntensity_baseline,numtrials,1);
surroundInten_baseline = repmat(surroundInten_baseline,numtrials,1);
maxInten_baseline = repmat(maxInten_baseline,numtrials,1);
width_baseline = repmat(width_baseline,numtrials,1);



vesselSegmentID = [1:numTrackedSegs]';vesselSegmentID = repmat(vesselSegmentID,numtrials,1);
stimX = uid(stimID_idx,1);
stimY = uid(stimID_idx,2);
stimClass = uid(stimID_idx,3);
stimSize = uid(stimID_idx,4);
stimX = repmat(stimX',numTrackedSegs,1);stimX = stimX(:);
stimY = repmat(stimY',numTrackedSegs,1);stimY = stimY(:);
stimClass = repmat(stimClass',numTrackedSegs,1);stimClass = stimClass(:);
stimSize = repmat(stimSize',numTrackedSegs,1);stimSize = stimSize(:);
temp = cell(size(totalgcamp));
temp(:) = cellstr(mouseID);
mouseID = temp;
sessionID = repmat(sessionID,size(totalgcamp,1),1);
trialnum = 1:numtrials;trialnum = repmat(trialnum,numTrackedSegs,1);trialnum = trialnum(:);
r = table(mouseID,sessionID,trialnum,stimX,stimY,stimClass,vesselID,prestimgcamp,totalgcamp,distance2gcamp,evokedMaxVessel,evokedIntegratedMaxVessel,evokedMaxSurround,...
    evokedIntegratedSurround,evokedMaxTotalVessel,evokedIntegratedTotalVessel,vesselSegmentID,stimSize,...
    vesselTtotalIntensity_baseline,surroundInten_baseline,maxInten_baseline,width_baseline,time2peak);

if dosave
    savename = strrep(trackingfn,'.mat','_resultsTable.mat');
    save(savename,'r');
end